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Artificial  Cooling  Due  to  Quiet  Injection  in 
Bounded  Plasma  Particle  Simulations 

William  S.  Lawton 
13  May  1987 

ABSTRACT 

An  explanation  is  proposed  for  an  artificial  cooling  effect  seen  in  electrostatic  particle- in-cell 

plasma  simulations.  Further  simulations  are  done  which  test  and  support  the  explanation. 

Introduction 

The  system  in  which  the  effect  was  observed  is  meant  to  simulate  a  Q- machine.  Electrons 
and  ions  are  injected  Cram  a  hot  plate  (the  left  side  of  the  system,  which  is  chosen  as  x  =  0), 
and  particles  which  strike  either  side  of  the  system  are  absorbed.  In  this  case,  the  end  plates  are 
electrically  isolated  (an  open  external  circuit).  The  density  of  injected  ions  is  chosen  to  be  much 
larger  than  the  density  of  injected  electrons,  so  that  a  potential  maximum  is  created  near  the  hot 
plate,  thus  creating  a  trapping  well  for  electrons.  The  simulation  parameters  are  shown  in  Table  1. 

Recent  simulations  of  Q-machines  [1]  have  shown  a  surprising  effect:  the  cooling  and  increase  in 
density  of  electrons  which  are  trapped  in  the  sheath  near  the  emitting  (hot)  plate.  Electrons  which 
are  trapped  in  this  potential  well  are  found  to  cool  over  long  time  scales  to  less  than  a  third  of  the 
temperature  at  which  they  were  emitted. 

Initially  the  simulation  region  is  empty  of  particles.  At  t  *  0,  ions  and  electrons  begin  entering 
the  system,  and  within  a  few  ion  transit  times  establish  a  near-equilibrium  which  agrees  with  what 
Vlasov  theory  suggests  the  equilibrium  state  should  be.  The  potential  profile  for  this  near-equilibrium 
(from  the  simulation)  is  shown  in  Fig.  1.  The  electron  and  ion  phase  speces  for  this  same  time  are 
shows  in  Fig.  2.  Next,  over  a  long  time  scale  (tens  of  ion  transit  times),  the  number  of  electrons 
in  the  trapping  well  steadily  increases,  the  temperature  of  the  trapped  electrons  decreases,  and  the 
spatial  width  of  the  trapping  well  increases  (see  Figs.  3  and  4).  This  process  has  not  saturated  in 
any  of  the  runs  which  will  be  shown  here,  but  saturation  does  occur  later,  at  a  trapped  electron 
temperature  which  is  much  smaller  than  the  injection  temperature  and  a  trapping  well  width  which 
spans  most  of  the  simulation  region. 
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No  simple  model  predicts  this  cooling.  The  Vlasov  model  (fully  time-dependent)  predicts  that 
the  phase  space  density  f(x,  v,  t)  in  the  trapped  well  may  be  no  higher  than  the  largest  value  of 
f(x,  v ,  t )  at  the  point  of  injection  (x  =  0  for  our  simulations).  For  the  Q- machine,  this  is  also  the 
value  of  /  on  the  separatrix  which  defines  the  trapping  region  of  phase  space.  Alternatively,  one 
might  expect  that  the  distribution  in  the  trapping  well  would  fill  in  through  some  scattering  process 
(due  to  the  discrete  particle  model)  to  a  Maxwellian  shape.  In  fact,  the  distribution  fills  in  far  past 
this  point. 

Some  possible  explanations  for  this  effect  are:  numerical  inaccuracy,  enhanced  fluctuations  due 
to  an  instability,  and  collisions.  None  of  these  explanations  seem  capable  of  accounting  for  the 
observed  effect.  The  consideration  of  fluctuations,  however,  led  to  what  seems  to  be  the  correct 
explanation. 


An  enhanced  level  of  fluctuations  would  be  expected  to  increase  particle  diffusion  in  velocity 
space,  but  in  such  a  way  as  to  heat  trapped  electrons.  In  thermodynamic  equilibrium,  this  heating 
effect  is  balanced  by  a  drag  on  the  particles,  giving  rise  to  the  Maxwellian  distribution.  The  level 


a f  fluctuations  in  thermodynamic  equilibrium  is  not  aero,  so  an  abnormally  low  level  of  fluctuations 
can  be  expected  to  cool  electrons.  Such  a  depressed  level  at  fluctuations  is  in  fact  specifically 
introduced  in  order  to  in  crease  the  signal-to-notse  level  of  most  simulations.  This  artificially  low 
level  of  fluctuations  appears  to  be  the  cause  of  the  artificial  cooling  observed  in  the  Q-machipe 
simulations.  In  thirty  years  of  periodic  rimnlations,  artificial  cooling  like  this  has  not  been 

reported.  A  reeson  for  this  effect  going  unnoticed  will  also  be  given. 

Ruling  Out  Some  Possible  Causes 

Three  candidates  far  the  cause  of  the  cooling  effect  can  be  eliminated.  These  are  numerical 
inaccuracy,  fluctuations,  and  collisions 

Numerical  inaccuracy  was  tested  by  <•*»*« ging  the  simulation  parameters  Ax  and  At,  which 
rsprawnt  the  x*>iml  find  (which  limits  the  spatial  resolution),  and  the  time  step  (which 


limits  the  temporal  resolution).  When  one  or  both  of  these  parameters  was  reduced  (improving  the 
resolution),  the  cooling  effect  remained  with  the  same  magnitude. 

An  enhanced  level  of  fluctuations  can  be  because  they  are  a  heating  influence.  In- 

trussing  the  fluctuation  level  enhances  diffusion  in  velocity  in  both  directions,  and  so  should  create 
just  ss  much  flux  out  of  the  trapping  well  as  into  it.  In  fact,  when  the  density  of  particles  in  phase 
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■pace  within  the  well  is  higher  that  outside  it,  one  would  expect  more  particles  to  leave  the  well 
than  enter  it. 

Collisions  are  harder  to  eliminate.  In  three  dimensions,  a  Fokker-Planck  analysis  (see  Ap¬ 
pendix  A)  shows  that  ions  moving  above  a  certain  z-velodty  (approximately  \/3  times  the  thermal 
speed  of  ions  of  the  same  temperature  as  the  electrons)  heat  the  electrons,  while  ions  moving  slower 
than  this  velocity  cool  them.  (The  Fokker-Planck  model  is  of  questionable  validity  here  due  to  the 
short  scale  length  of  the  trapping  well,  but  it  should  be  indicative.)  Thus,  while  the  ions  between 
the  hot  plate  and  the  center  of  the  well  should  do  some  cooling,  the  ions  an  the  other  side  of  the 
center  of  the  well  should  do  just  about  as  much  heating.  The  electrons  which  are  not  trapped  should 
also  do  some  coding  throughout  the  well,  but  the  effective  temperature  of  the  passing  electrons  is 
very  to  the  emission  temperature,  so  the  strong  coding  of  the  trapped  electrons  cannot  be 
explained  by  electron-electron  rdlisions.  The  overall  balance  of  these  effects  is  difficult  to  predict, 
bat  is  likely  to  be  small,  and  the  equilibrium  state  would  dearly  be  much  closer  to  the  collisionleas 
state  than  is  observed  in  simulation. 

In  one  dimension,  however,  the  rolliaton  pro  cam  is  very  different  [2].  For  instance,  particles  of 
the  same  aperies  do  not  scatter  when  they  collide,  they  either  pans  through  each  other,  or  exchange 
velocities.  Thus,  the  passing  elsctrons  cannot  be  cooling  the  trapped  electrons  through  cdliaioos 
in  this  simulation.  Part  ides  of  differing  species  do  not  collide  in  the  usual  sense  of  introducing  a 
random  change  in  velocity,  either.  Their  mBirione  exchange  energy  in  n  very  rimple  and  much  more 
predictable  way  —  either  fay  reflecting  velocities  in  the  center  of  mane  frame,  or  paming  through  each 
other  with  unaltered  velocities.  This  has  the  remit  of  greatly  diminishing  the  effects  of  collisions 
relative  to  the  three  dimension  si  case.  Given  these  facts  in  combination  it  seems  implausible  that 
coOirioua  could  be  responsible  for  the  effect.  This  is  admittedly  a  weak  argument;  the  strongest 
evidence  for  nrllieinne  being  unimportant  is  the  rimnlatioo  evidence  presented  here  for  the  depressed 
level  of  fluctuations  being  the  only  pro  ram  of  importance. 

Electric  Arid  fluctuations 

Thermal  energy  may  reride  in  one  of  two  reservoirs  in  a  plasma-  disorganised  particle  motion 
and  fluctuations  of  the  electrostatic  wave  field.  The  word  field  is  used  here  not  in  the  sense  of  the 
electric  field,  but  in  the  sense  of  the  set  of  waves  which  are  supported  by  a  medium  (as  in  Quantum 
Field  Theory).  It  m  necessary  to  consider  the  electrostatic  wave  field  as  being  something  sepsrste 


from  the  electric  field,  since  the  organized  motion  of  particles  plays  a  role  in  the  electrostatic  waves  in 
the  plasma  At  short  wavelengths  the  distinction  between  wave  fluctuations  and  particle  fluctuations 
becomes  Murry,  but  at  long  wavelengths  the  distinction  is  clear. 

There  is  coupling  between  these  reservoirs,  and  energy  flows  between  the  two.  In  thermodynamic 
equilibrium,  the  flow  of  energy  from  wave  fluctuations  to  particles  is  exactly  balanced  by  the  flow  of 
energy  from  particles  to  wave  fluctuations.  These  flows  are  called  respectively,  Landau  damping  and 
Cerenkov  emission,  and  the  flow  rates  are  proportional  to  the  energy  present  in  the  reservoir  which 
is  losing  energy  (see  Fig.  5).  The  balance  between  these  flows  in  thermodynamic  equilibrium  gives 
rise  to  the  Fluctuation-Dissipation  Theorem,  which  relates  the  temperature  to  the  level  of  electric 
field  fluctuations  in  the  p1 

Quint  injection 

In  order  to  reduce  the  background  noise  level  in  particle  simulations,  particles  are  usually  put 
into  phase  space  (either  loaded  initially  or  injected  over  time)  nearly  uniformly.  This  has  the  effect  of 
drastically  reducing  the  amplitude  of  fluctuations  in  the  electric  field  [3].  For  computer  runs  of  short 
duration,  this  is  entirely  beneficial.  Even  for  long  runs,  if  the  boundary  conditions  are  periodic,  the 
worst  that  can  happen  is  that  the  fluctuation  level  of  the  electric  Arid  ((£*))  rises  up  to  the  natural 
lead,  and  since  the  total  boat  content  of  the  fluctuations  is  not  large  (».«.,  ta(E?)/nkT  is  small),  and 
the  system  is  dosed,  the  particles  do  not  cool  appreciably  (see  fig.  6).  This  cooling  of  particles  in 
periodic  simulations  has  not  been  explicitly  reported,  but  should  be  observable.  (Simulation  experts 
may  be  reminded  of  Gitomar  and  Adam  [4],  but  the  effect  they  noted  was  a  rapid  rise  in  the  level 
of  fluctuations  doe  to  a  multibeam  instability  —  a  very  different  phenomenon  from  that  described 
here.) 

A  bounded  aimn latino  is  an  open  system,  and  the  situation  is  quite  different.  Most  particles 
are  in  the  system  only  for  one  or  two  transit  times.  Only  the  particles  which  are  trapped  stay  for 
long  times.  The  wave  fluctuations,  however,  are  not  trapped  at  all.  Since  they  are  collective  inodes, 
they  move  at  their  group  velocities  (which  are  strongly  affected  by  the  bulk  mot  too  of  the  plasma). 
Thus  the  wave  fluctuations  can  and  do  leave  the  system,  only  to  be  replaced  with  the  artificially 
reduced  fluctuations  in  the  newly  injected  particles.  In  other  words,  the  trapped  electrons  can  cool 
via  Cerenkov  rmissinn  of  waves  which  then  leave  the  system.  The  energy  flow  is  now  much  more 
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complicated  than  in  the  periodic  caae  (see  Fig.  7).  This  mechanism  is  completely  physical;  if  the 
noise  level  of  injected  electrons  could  be  reduced  in  the  laboratory,  this  effect  would  occur. 

Simulation  evidence 

Simulation  offers  many  methods  of  testing  the  hypothesis  that  the  quiet  injection  of  particles 
is  responsible  for  the  cooling  of  the  trapped  electrons.  Several  of  these  were  tried,  and  all  supported 
the  hypothesis. 

Figure  8  shows  the  results  of  a  quiet  start  run.  Plotted  are  the  electron  velocity  distribution 
function  averaged  over  a  narrow  region  near  the  potential  maximum  at  the  end  of  the  run,  and  the 
time  history  of  the  total  number  of  electrons  in  the  system.  Note  the  elevated  electron  density  and 
the  decreased  temperature  of  the  trapped  electrons  in  the  distribution  function  plot,  and  the  steady 
increase  in  the  total  number  of  electrons. 

The  most  obvious  test  of  the  hypothesis  is  to  replace  the  quiet  injection  scheme  with  one  which 
has  full  noise,  meaning  that  the  times  of  injection  are  random  with  a  uniform  distribution  in  time, 
and  that  the  velocities  of  injection  are  random  with  a  half-Maxwellian  distribution.  The  result, 
shown  in  Fig.  9,  is  that  the  system  settles  down  to  an  equilibrium  very  quickly,  and  the  trapped 
electrons  never  cool  below  the  temperature  of  injection. 

One  more  piece  of  information  can  be  inferred  from  this  run.  A  valid  objection  to  the  support 
this  run  gives  the  hypothesis  is  that  any  source  of  noise  will  cause  heating,  and  so  it  is  not  obvious 
that  the  correct  source  of  noise  has  been  found  or  that  lack  of  noise  is  indeed  responsible  for  the 
effect-  This  run  shows,  however,  that  the  noise  is  of  roughly  the  right  magnitude.  If  the  noise  level 
were  too  high,  then  the  trapped  electrons  in  the  noisy  case  would  reach  a  Maxwellian  profile  sooner 
than  the  quiet  case  (which  initially  comes  to  the  flat  distribution  predicted  by  Vlasov  theory),  since 
fluctuations  enhance  diffusion  in  both  directions  across  the  separatrix.  If  the  noise  level  were  too 
low,  some  residual  moling  effect  would  be  seen.  That  neither  of  these  effects  are  seen  lends  support 
for  the  hypothesis. 

To  further  test  the  hypothesis,  some  runs  were  made  with  four  tim—  as  many  particles.  Since 
the  physical  parameters  were  kept  constant,  the  charge  and  mass  of  each  particle  was  decreased  by  a 
factor  of  four,  la  the  quiet  start  case,  since  the  electrons  are  injected  quasi- regularly,  the  fluctuations 
in  the  charge  density  should  vary  inversely  with  the  number  of  particles,  in  this  case  a  factor  of  four. 
This  in  turn  should  result  in  a  factor  of  four  reduction  in  the  RMS  fluctuation  of  the  electric  field 
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(( E  —  E)3)1^3.  The  diffusion  rate  is  proportional  to  {(E  —  E)3),  so  the  rate  at  which  each  particle 
diffuses  into  the  trapping  well  should  be  reduced  by  a  factor  of  16.  (Since  there  are  four  times  as 
many  particles,  the  absolute  number  of  particles  diffusing  into  the  trapping  well  should  be  reduced 
by  a  factor  of  four.) 

In  the  case  of  full  thermodynamic  noise,  the  Fluctuation-Dissipation  Theorem  dictates  that 
increasing  the  number  of  particles  by  a  factor  of  four  must  result  in  the  mean-square  electric  field 
fluctuation  {(E — E)3)  decreasing  by  a  factor  of  four.  This  in  turn  should  reduce  the  rate  of  diffusion 
by  only  a  factor  of  four.  (In  this  case,  the  absolute  number  of  particles  diffusing  into  the  trapping 
well  in  a  given  time  interval  will  not  be  reduced.)  As  Figs.  10  and  11  show,  the  quiet  and  noisy 
runs  obey  roughly  the  expected  behavior.  The  electric  field  was  also  followed  to  see  if  it  obeys  the 
behavior  described  in  the  last  paragraph,  and  as  Fig.  12  shows,  it  does. 

Conclusion 

Simulations  support  the  hypothesis  that  quiet  loading  in  bounded  plasma  particle  simulations 
causes  artificial  cooling  of  trapped  electrons.  This  effect  is  very  small,  and  therefore  has  not  been 
reported,  in  periodic  simulations  because  of  the  low  heat  capacity  of  the  wave  fluctuation  field,  but 
it  will  be  a  limitation  on  the  accuracy  of  some  bounded  plasma  simulations. 
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Appendix  A:  Fokker-P lanck  Analysis  of  Heating 


To  find  the  heating  effect  of  ions  on  the  trapped  electron  distribution  in  the  three-dimensional 
case,  the  Fokker-Planck  equation  can  be  used.  I  will  consider  a  spherically  symmetric  electron 
distribution  (later  a  Maxwellian),  and  a  cold  beam  of  ions  with  which  the  electrons  will  collide.  (The 
solution  to  this  problem  can  be  used  as  a  Green's  function  for  more  complicated  ion  distributions.) 
Since  the  electrons  all  by  themselves  should  move  toward  a  Maxwellian,  the  effect  of  collisions 
between  electrons  should  not  be  of  interest  in  the  context  of  cooling.  Naturally,  with  a  beam  of  ions 
moving  through  it,  a  spherically  symmetric  distribution  will  not  stay  spherically  symmetric  for  long, 
but  the  collisions  among  electrons,  and  the  trapping  well  in  which  the  electrons  reside  will  tend  to 
return  the  distribution  function  to  being  nearly  spherical.  Also,  we  are  not  interested  in  a  complete 
analysis,  but  only  whether  the  beam  imparts  energy  to  the  electrons  or  takes  it  from  them,  i.e., 
whether  dT/dt  is  positive  or  negative. 

The  Fokker-Planck  equation  can  be  written  [5] 

*  -r*e  •  {/V.h  -\v9.  (/V.V.9)  j  (al) 

where 

Vah  *  -*»“//  *  -  Vo)  («*2) 

and 

V9gm?th  (o3) 

m 

with  T  being  a  constant  whose  magnitude  is  unimportant  here,  //  representing  the  distribution 
function  of  the  field  particles  (in  this  case  a  cold  beam  of  ions),  and  ft  representing  the  electron-ion 
reduced  mass, 

mM 

The  equations  for  h  and  g  can  be  solved  to  yield 


and 


h 


y  ■  (r  -  uj| 


but  these 


will  not  be  needed  in  the  following  calculation. 
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FYom  (al),  the  time  derivative  of  the  kinetic  energy  due  to  collisions  with  ions  is 

^  |  5  =  -I'm  J  ^V.  •  j/V.h  -  iv.  •  (/V.V.p)}  (a4) 

Integrating  the  right  hand  side  by  parts  (once  for  the  first  term  and  twice  for  the  secon  term) 
gives 

^  »  Tm  J  {«/•  fV.h  +  i/Vjp}  d»u 

=  Tfif  •  («/)  +  /J  hfv  (o5) 

Now  writing  m/p  as  1  +  m/M, 

^  =  -r^y  (£v.  •  (vf)  +  (2/ + «f  •  V./)}  /.<*»  t,  (06) 

It  is  now  possible  to  avoid  a  lot  of  algebra  by  defining 

m-f  »'/(v')dw'  (o7) 

and 

v>2f(v')dv'  («*) 

v  Jo 

(this  is  where  the  assumption  of  spherical  symmetry  is  used).  We  then  note  that 

V.p  =  -vf  (a9) 


When  1*0  is  small,  p  dominates  since  it  goes  to  a  positive-definite  constant  whereas  q  falls  off 
as  Vq.  When  1*0  is  large,  q  dominates,  since  it  falls  off  as  1/vq  and  p  falls  off  like  /  —  typically 
exponentially.  Thus  the  expected  behavior  is  reproduced;  slow  particles  cool  a  distribution,  and  fast 
ones  heat  it.  The  transition  point  is  the  solution  of 

?(« o)  =  ~£p(v o)  («13) 


Now  let  us  confine  ourselves  to  a  Maxwellian  distribution  to  find  an  approximation  to  the 
velocity  which  neither  beats  nor  cools  the  distribution.  Assuming  m/M  is  reasonably  small,  vo 
should  be  small,  so  p  and  q  can  be  expanded  in  a  Taylor  series.  For  a  Maxwellian  distribution, 
p  =  v*f,  so  for  small  v,  p  — ►  vJ/( 0)  to  leading  order.  For  any  distribution  at  all, 


Setting  these  equal  gives 


*§  *  =  3vf 


(*14) 


This  approximation  is  good  to  one  percent  even  for  mass  ratios  as  high  as  1/40. 
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Fig.  2.  Electron  (a)  and  ion  (b)  phase  spaces  for  potential  in  Fig.  1. 
Note  smoothness  and  patterns  due  to  quiet  injection. 
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Fig  5  Energy  balance  between  particle  thermal  energy  and  wave  fluctuations  in  periodic  or  infinite  plasma 
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Fig.  6.  Transition  from  initial  quiet  start  (a)  to  final  equilibrium  (b)  takes  place  on  a  slow 
time  scale  The  final  energy  content  of  the  waves  is  much  smaller  than  that  of  the  panicles. 
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Fig.  12  Time  histones  of  etectnc  held  in  center  of  simulation  region  for  four  different  runs 
(ai  and  (b>  art  quiet  injection  runs  with  <b)  having  four  times  as  many  particles  as  (a),  and  tc v 
and  (d)  are  noisy  injection  runs  with  (d)  having  four  times  as  many  panicles  as  (c)  RMS  value 
of  £  in  (bl  should  be  1/4  that  of  (a),  and  RMS  value  of  £  in  (d)  should  be  half  that  of  (c) 


